The present invention relates to a method for estimating aspects of the bone quality or skeletal status of a vertebrate on the basis of two-dimensional image data comprising information relating to the trabecular structure of at least a part of a bone, typically a long or tubular bone. Especially bones comprising a tube (diaphysis) and at least at one end an epiphysis comprising trabecular bone are suitable for use in the present method.
In the present context, skeletal status is taken to relate to e.g. biomechanical competence, trabecular homogeneity, cortical thickness of the bone or a quantification of sensitivity, specificity or accuracy of diagnosis of osteoporosis.
Age, or more likely co-variates with age, such as e.g. risk of falling, the type of falls and skeletal status, significantly affects fracture risk. Hence, any type of BMD must be viewed only as an index of risk. However, when age-matched individuals with and without fractures are compared, an overlap in BMD measurements is apparent. It is contemplated that this, in part, is due to the fact that important skeletal properties, such as trabecular connectivity (or homogeneity), is largely ignored by BMD measurements, since small losses in cancellous bone mass may lead to a disproportinately large loss in biomechanical competence.
The surface area of cancellous bone is greater than that of cortical bone, even though cancellous bone may occupy only 25% of the skeletal mass of a healthy individual and as little as 10% of that of an osteoporotic patient. Because of the high surface-to-volume ratio of cancellous bone tissue, osteoporosis is thought to affect trabecular sites early in the disease process. In postmenopausal osteoporosis, there is a marked loss of trabecular elements, but little evidence of actual trabecular thinning. The loss of elements is thought to be a result of cavities that transect or perforate structures. This is thought to have important implications for skeletal strength, and hence for those components of fracture risk related to skeletal strength.
Currently, the connectivity of cancellous bone can be measured reliably only by invasive techniques or by the technique described in the Applicants co-pending International application No. PCT/DK95/00338.
In U.S. Pat. No. 4,903,203, a method of positioning an ROI in a bone is disclosed. In this reference, only a single ROI is determined, mined, and no variance of PTD is mentioned.
In EP-A-0 648 467, different methods are disclosed wherein BMD is evaluated in a distal radius. In order to evaluate the method herein, BMD was investigated at different positions along the length of the bone. However, the actual method of this reference relates to a single BMD value averaged over a number of measurements. Thus, no variation of PTD has been investigated or is pointed to.
Presently, a new method is proposed that by quantifying trabecular homogeneity, as it appears on e.g. a radiograph of the distal forearm, an index of skeletal status may be provided.
Researchers have proposed a number of different ways to measure variations in radiographs that can be attributed to cancellous bone. Geraets et. al. (W. G. Geraets, P. F. Van der Stelt, C. J. Netelenbos and P. J. Elders, xe2x80x9cA new method for automatic recognition of the radiographic trabecular patternxe2x80x9d, J. Bone Miner. Res. 5(3), pp 227-233, March, 1990) proposed different textural measures estimated from a 10 mmxc3x9710 mm Region Of Interest (ROI), extracted close to the border of the former epiphysal growth plate. The extracted measures were correlated to wrist BMD, and moderate correlations were obtained. The inventor has not been able to obtain statistically significant correlations with biomechanical competence (of the hip), using this method. Weinstein et. al. (R. S. Weinstein and S. Majumdar, xe2x80x9cFractal geometry and vertebral compression fracturesxe2x80x9d, J. Bone Miner.. Res., 9(11), Nov. 1994) and, Majumdar et. al. (S. Majumdar, R. S. Weinstein and R. R. Prasad, xe2x80x9cApplication of fractal geometry techniques to the study of trabecularbonexe2x80x9d, Med. Phys., 20(6), pp 1611-1619, November 1993) proposed applying a box counting algorithm for fractal analysis of the trabecular bone matrix.
Chung (H. W. Chung. C. C. Chung, M. Underweiser and F. W. Wehrli, xe2x80x9cOn the fractal nature of trabecular structurexe2x80x9d, Med. Phys., 21(10), pp 1535-1540, October 1994) investigated this method and concluded that the box counting algorithm is not suitable for fractal analysis of cancellous bone and continues to suggest that the fractal appearance of the trabecular network is artifactual. Benhamou et al. (C. L. Benhamou et al., xe2x80x9cFractal organization of trabecular bone images on calcaneous radiographsxe2x80x9d, J. Bone Miner. Res., 9(12), pp 1909-1918, December 1994) suggests that the surface may, in fact, be fractal, if it is not binarized. It is proposed to use a fractional gaussian noise model to obtain an estimate of the fractal dimension for a number of calcaneous images, and a high degree of reproducibility of this method is reported, although correlations to specific biomechanical properties are not reported in this paper.
Trippi et al. (D. Trippi, M. Chimenti and R. Bozzi, xe2x80x9cA computer-assisted method for the study of the trabecular . . . xe2x80x9d, J. Digit. Imaging, 6(2), pp 140-147, May 1993) propose measuring a coarseness index, based on image variance, on the forearm, and find a decrease in coarseness in the distal-proximal direction. The paper aims at defining a systematic way of selecting an ROI, and recommends measuring 6 to 12 mm from the former epiphysal growth plate.
According to the invention, it has been found that it is of limited use to look at one particular Region Of Interest. Instead, changes in the local variations, attributable to cancellous bone, must be quantified in order to obtain meaningful information about trabecular homogeneity. These findings also suggest that very strong correlations between BMD and measures of trabecular homogeneity may not exist, due to very fundamental differences in the nature of such measures. Finally, trabecular homogeneity measured on eg the forearm offers important independent information about biomechanical competence of the hip, that is not offered by any densitometric measures obtained on the appendicular skeleton.
Consequently, the present invention relates to a method for estimating the skeletal status of a vertebrate on the basis of two-dimensional image data comprising information relating to the trabecular structure of at least a part of a bone of the vertebrate, the image data being data obtained by exposing at least the part of the bone to electromagnetic radiation, the method comprising subjecting the image data to a statistical analysis comprising:
a background correction procedure in which low frequency intensity variations not related to the trabecular structure of the bone is reduced relative to image data related to the trabecular structure of the part of the bone,
a feature extraction procedure comprising
(a) determining values reflecting the projected trabecular density in the image data, caused by the X-ray attenuating properties of cancellous bone in the part of the bone, for each of a number of locations or areas in the image data,
(b) deriving one or more features from the variation of the determined PTD-values, preferably in the longitudinal direction of the bone, and
an estimation procedure in which the skeletal status of the vertebrate is estimated on the basis of the one or more derived features and optionally other features related to the bone of the vertebrate and a predetermined relationship between the features and reference skeletal status parameters.
As indicated above, in the present context xe2x80x9cskeletal statusxe2x80x9d is not equalled to xe2x80x9cbone quantityxe2x80x9d, such as Bone Mineral Density, as a loss of trabecular elements may lead to a disproportionate loss of biomechanical competence.
In fact, it is contemplated that a better estimate of skeletal status may be obtained by considering pertinent properties of cortical and cancellous bone individually instead of treating these two types of bones similarly, as it is done in standard densitometric methods.
The number of areas or positions may be any number higher than one. Normally, the fewer areas or positions, the more difficult will it be to provide a robust and reproducible method, as eg the extent to which the cancellous bone of a healthy individual extends away from the distal end of eg a distal radius may vary quite a lot. Therefore, it is preferred that the number of areas or positions is larger than 10, such as larger than 100, preferably larger than 500.
An X-ray image is a projection of the attenuating properties of the tissue in the direction of the X-ray radiation and is thus a projection of the 3-dimensional distribution of the X-ray attenuating properties of the cancellous tissue. The density of this projected structure yields a projected trabecular density. Thus, even though the bones and the trabecular structures thereof are inherently three-dimensional, the projection of this structure onto two dimensions, such as in a radiographic image, conveys useful information about the spatial distribution of cancellous bone.
It may be preferred to not only determine a single value for each location or area corresponding to local variation. In fact, two or more values reflecting different aspects of the projected trabecular density may be determined for each location or area.
Determining, evaluating, predicting or using the PTD at the locations or areas of the bone will, naturally, provide information relating to the trabecular homogeneity of the bone. As is contemplated in the present invention, this trabecular homogeneity reflects the skeletal status.
Low frequency intensity variations may be due to, e.g., scattering of the electromagnetic radiation, anatomic structures surrounding the illuminated part of the trabeculae such as cortical bone, fat tissue and muscles of varying thickness as well as, e.g., inhomogeneous X-ray illumination of the part of the bone. As these features are irrelevant to quantification of trabecular homogeneity, the low frequency intensity variations are preferably removed from the image data.
The features extracted from variations in the PTD-values are used for calibrating a model that describes any given reference variable relating to the skeletal status as previously described. Subsequently, the calibrated model may be used for predicting pertinent expressions of skeletal status for a given vertebrate.
As described above, bone quality or skeletal status is not equal to bone quantity in the present context, and the reference skeletal status parameters are preferably parameters related to the biomechanical competence of the bone, sensitivity, specificity or monitoring of treatment effects rather than, e.g., BMD or BMC.
However, to provide even higher diagnostic sensitivity, combination of the present invention with a measurement of BMC/BMD is contemplated to be useful. BMC/BMD measurements an be obtained either by the well-known single or dual wavelength absorptiometry technique or they may, if the image data are obtained on the basis of an X-ray image, be available from the same image data, where, e.g., a standard aluminum wedge has been illuminated together with the bone in question.
At present, the most common way of obtaining non-invasive information about the bone quantity and skeletal status is based on densitometric methods. According to the invention, the image data are preferably obtained for X-ray images. X-ray images may easily be used to generate the image data for use in the present invention. In order to be able to use image processing, an essentially analogue X-ray image is preferably digitized and introduced into a computer by scanning the X-ray image. Alternatively, methods using direct measurement of the X-rays may be used.
Naturally, it is preferred that the scanning of X-ray images is performed with a sufficiently large resolution and dynamic range in order to ensure that a minimum relevant information of the X-ray image is lost in the transfer to the image data. Thus, it is preferred that the scanning has been performed at a resolution of at least 4 pairs of lines per centimetre, such as at least 5 pairs of lines per centimetre. Even though this may be sufficient for obtaining a satisfactory estimation of the skeletal status, it is preferred that the resolution of the scanning of the X-ray film is at least 10 pairs of lines per centimetre, such as at least 25 pairs of lines per centimetre, preferably at least 50 pairs of lines per centimetre. It is possible to obtain resolutions in the scanning up to at least 100 pairs of lines per centimetre, such as at least 250 pairs of lines per centimetre and probably as high as 500 pairs of lines per centimetre, such as at least 600 pairs of lines per centimetre.
Furthermore, it is preferred that the resolution of the scanner is better than 4 true bits per pixel, such as better than 6 true bits per pixel, preferably better than 8 bits per pixel, such as better than 12 true bits per pixel, preferably better than 16 bits per pixel.
Naturally, in order to have the scanner actually scan the image, it should be able to transilluminate the radiographs.
The ESKOSCAN 2450 from Eskofot A/S has been found to fulfil the above criterias and to be highly suited for use in the method of the invention.
According to the invention, evaluating the homogeneity of the trabecular bone is performed by determining the local variation in density caused by variations in the cancellous bone composition resulting in a variation in the projected trabecular density seen in the image data at a number of locations or areas in the image data, preferably in a longitudinal direction of the bone. Subsequently, one or more features for use in the estimation procedure are derived from the PTD values.
Even though it is contemplated that the variation of PTD may be used in the present method for substantially all directions in the bone, especially that of the longitudinal direction is found interesting, and especially in long bones, such as the tibia, ulna, and the distal radius, where strong relations to other measures of skeletal status have been found.
In fact, it is contemplated that a large number of features derived from this variation of PTD in the longitudinal direction of the bone will provide a significant correlation with the failure load of e.g. the hip.
In order to facilitate the determination of the Regions Of Interest used in the bones and to enhance repeatability of the method, it is preferred that the PTD-values are determined at locations or areas situated longitudinally spaced apart within the bone. It is currently preferred that the locations or areas for which the PTD-values are determined are located at predetermined areas or positions in the image, such as along a reproducibly positionable curve in the image, preferably a substantially straight line substantially along a longitudinal axis of the bone.
The PTD-values may be determined in a number of ways, and they may be chosen to relate to different properties of the projected trabecular structure. At present, the PTD-value determined for each of the locations or areas is a value relating to first, second, third or higher order properties of the image data relating to the projected trabecular structure at the location or in the area in question.
Preferably, the PTD-values are calculated on the basis of a determination of any property of the Fourier power spectrum, a parametric spectral estimate or of the grey level co-occurrence matrix of information in the image data at the location or in the area in question, the variance of the image data at the location or in the area in question, or a mean value of the size and density of edges in the image data at the location or in the area in question.
Preferably, a PTD-profile is defined by the determined PTD values and a (preferably reproducible) measure relating to the positions of the areas or positions in relation to each other or in relation to a fixed position within the bone, such as, the determined PTD-values as a function of eg a distance from the location or area for which the individual PTD-value was determined to a predetermined reference point. In this situation, at least one of the above one or more features is derived from the profile. Features of this profile have been found to correlate especially well with parameters of skeletal status. Preferably, the predetermined reference point is also located on the line along which the PTD values are extracted, and preferably below the former epiphysal growth plate.
Thus, the PTD-profile preferably describes the PTD as a function of the distance between the area or position at which the individual value was derived and the predetermined reference point.
This profile has the advantage that reproducible information may be derived from all bones independently of the size and shape of the bone in question. Moreover, this type of measurement gives a global description of the skeletal status as opposed to measuring in a region of interest chosen more or less ad hoc.
If, as is preferred, the locations or areas for which the PTD-values are determined are situated between a former epiphysal growth plate of the bone and a second location within the bone at which the image data show substantially no projected trabecular structure, the profile will reflect a falling PTD, as the distance from the epiphysal growth plate increases.
It is of interest to consider a maximum and/or a minimum of the profile or a smoothed version thereof. Features may therefore also be derived from these measures.
The fluctuation or variation of the PTD-profile is also of interest. Therefore, one or more features are preferably derived from this fluctuation or variation. Such a feature may be derived from a polynomial fitted to the profile or from the standard deviation of the difference between the polynomial and the profile or a smoothed version thereof. In these situations, the polynomial is preferably of an order in the interval 1-10, such as 2-5, preferably 3-4.
In addition, a feature may be derived from a steepness of the profile, as this steepness may describe the decrease in trabecular bone from the epiphysal growth plate towards the center or middle of the bone.
An optional or additional feature may be derived from the standard deviation of a difference between the profile or a smoothed version thereof and a morphological opening or closing of the profile or of a smoothed version of the morphological phological opening or closing of the profile.
It is presently preferred that the background correction procedure comprises at least reducing or optionally removing low frequency information having a frequency significantly lower than the average spacing of the projected trabeculae. As described above, this low frequency content of the image data is typically caused by varying thickness of cortical bone, fat tissue and muscles of varying thickness as well as, e.g., inhomogeneous X-ray illumination of the part of the bone. Dual X-ray absorptiometry was invented to circumvent such problems relating to compensating for soft tissue. Also, other non-invasive image acquisition techniques, such as, e.g., MR and CT imaging may generate image data not automatically or not to the same degree xe2x80x9csuffering fromxe2x80x9d this type of undesired effect.
Even though it is difficult to exactly quantify the limit between undesired low frequency information and the desired higher frequency relevant information independently of the image resolution, it is presently preferred that information having frequencies half or less than the average spacing of the projected trabeculae is at least reduced or optionally removed. More preferably, information having frequencies being a quarter or less, such as a tenth or less than the average spacing of the projected trabeculae is at least reduced or optionally removed.
One preferred way of reducing or preferably removing the low frequency information is using a background correction procedure comprising generating secondary image data as a result of performing a median filtering with a predetermined kernel size and subtracting this result from the original image data. One of the advantages of using a median filter is that this operation is edge preserving. This technique is known as xe2x80x9cunsharp filteringxe2x80x9d.
Another way of reducing or preferably removing the low frequency information is using a background correction procedure comprising generating secondary image data as a result of performing a mean filtering with a predetermined kernel size and subtracting this result from the original image data. A mean filtering is typically much faster than the median filtering. However, the mean filtering is not edge preserving. Hence, the unsharp filtering may remove data we wish to include in the measurement.
It is typically preferred that the kernel size is at the most xc2xd of the image data, such as at the most xc2xc of the image data, preferably at the most {fraction (1/10)} of the image data, more preferably at the most {fraction (1/20)} of the image data.
A third way of reducing or preferably removing the low frequency information is using a background correction procedure comprising globally fitting a two-dimensional polynomial to the image data and generating a background corrected image on the basis of the residuals of the fitting procedure. Apart from potential difficulties in determining the optimal order of the polynomial, which the person skilled in the art will know, this method potentially offers an extremely fast background correction of image data.
It is contemplated that a suitable order of the polynomial may be at the most 15, such as at the most 10, more preferably at the most 5.
Preferably, subsequent so the background correction procedure, the grey level distribution of the image data is standardized by performing a histogram stretch in order for the distribution to match a Gaussian. This step has the advantage that eg grey levels of different images may now be compared directly. This step is standard in image processing of this type.
Another method of background correcting the image data is one wherein the total absorption of the cortical bone in the bone is estimated at the different areas or locations and subtracted from the image data.
As bones, at the positions thereof comprising trabeculae, also comprise cortical bone, the total absorption thereof is a combination of the absorption of these two components. The below method has been devised for separating the two contributions. It is clear to the skilled person, that information relating to the cortical absorption in the image data will typically have a frequency lower than the information relating to the trabeculae.
Consequently, in a fourth embodiment, the background correction comprises:
(I) assessing, for each location or area, a value relating to the X-ray absorption caused by cortical bone, and
(II) subtracting, for each location or area, the value of assessment (I) from the two-dimensional image data.
Preferably, this correction is performed on the raw image data as taken directly from a scanned X-ray image or directly from eg an X-ray detecting CCD.
Preferably, the assessment (I) is performed by:
(III) assessing the cortical thickness at one or more locations of one or more bones in the image data,
(IV) assessing the absorption caused by the assessed thickness of cortical bone in the image data at the one or more locations,
(V) for each area or position, assessing the thickness of cortical bone causing X-ray absorption, and
(VI) assessing, for each location or area, a value relating to the X-ray absorption caused by cortical bone on the basis of the assessments (III), (IV), and (V).
Preferably, the absorption caused by the cortical bone is expressed as a mean value of grey levels at the position in question.
If eg an X-ray image is taken of a hand comprising both number of metacarpal bones and part of the distal radius, the Combined Cortical Thickness may be determined for a number of the metacarpal bones as well as at positions in the distal radius, where no trabeculae are present.
From these CCT measurements and the assumption that the cortical bone has the same thickness around the circumference of the bone at the longitudinal position thereof a relation may be determined which relates the absorption of the cortical bone to the thickness thereof. The respective absorption information may be derived from the image at the position where the X-rays have experienced the absorption of the CCTxe2x80x94that is at the same longitudinal position in the bonexe2x80x94in the middle of the bone.
This method has the advantage that, now, absorption data may be derived about the cortical bone (and on the basis thereof about the trabecular bone) which is independent of the manner in which the image was taken, developed, illuminated, scanned, etcxe2x80x94as long these properties are at least substantially the same over the surface of the scanned image.
This data may now be transferred to the part of the distal radius comprising trabeculae. Assuming again that the thickness of the cortical bone is substantially constant around the circumference of the bone (the distal radius), the thickness of the cortical bone as determined from the image data (in the same manner as used at the determination of the CCT) may be used in a determination of the absorption in the image relating to cortical bone and to trabecular bone.
Having determined the thickness and, using the above-determined graph or relation, the absorption of the cortical bone, the xe2x80x9cdeviationxe2x80x9d therefrom may be assumed to relate to the absorption from the trabeculae. This absorption may be determined as a function of the position along the longitudinal position along the length of the bone or simply at predetermined positions or areas.
This calculated total absorption (or a related density predicted therefrom) may subsequently be used for deriving features for the quantitative determination of the skeletal status.
In addition to the above-mentioned features derivable from the image data, it may be preferred to input additional data relating to the vertebrate in question but which may not be derivable from the image data. This type of data may be information relating to the age and/or sex, and/or species, and/or race and/or the specific bone considered in the vertebrate, and/or a estimated Bone Mineral Density, and/or a estimated Bone Mineral Content.
Even though the DXA- or SXA-BMD determined using other equipment may be introduced in the estimation procedure, this measure may also optionally be determined from the image data. BMD is estimated by including data from a reference object in the exposure of the bone to the electromagnetic radiation and on the basis of the absorption of the electromagnetic radiation of the bone and of the reference object.
An especially preferred additional information to add to the features derived from the PTD-values is the thickness of the cortical bone of the same or another tubular bone of the vertebrate. Improved results are obtained, when indices relating to the thickness of the cortical shell of a bone, such as the Combined Cortical Thickness or a relative measurement as the Metacarpal Index, are used in the estimation procedure. This is due to the fact that the cortical and the trabecular bone tissue are thought to be metabolically different skeletal components which both contribute to the overall skeletal status. Quantifying both components will therefore most probably yield a more precise picture of the skeletal status.
The above-mentioned model may be determined in a number of ways. However, establishing a model of this type requires obtaining corresponding values of all relevant features and an end point variable relating to eg the biomechanical competence, the diagnostic sensitivity, specificity or accuracy in the diagnosis of osteoporosis which is to be estimated by the estimation procedure.
The system will require a calibration against suitable index of skeletal status. Examples of interesting skeletal status parameters are the biomechanical competence, and/or a Bone Mineral Density measurement, and/or a Bone Mineral Content measurement, and/or treatment effect and/or fracture incidence and/or a score value by a skilled radiologist. The method of the present invention will output a calibrated prediction of the chosen end point variable given a set of input variables: if the method is calibrated towards a strength parameter, the output of the method will relate to the strength of the bone in question.
It is expected that the present invention will be valuable in predicting the biomechanical competence through trabecular homogeneity and/or cortical thickness of the bone and/or improved sensitivity, specificity or accuracy in the diagnosis of osteoporosis.
As the calibration of the estimation procedure will reflect a predetermined relation between image features and skeletal status, this model should be chosen in a manner so that it correlates to the skeletal status. The assessment of the extent to which the embodiments of the present invention are able to estimate skeletal status may be performed by:
1. Prospective cohort-studies, in which odds ratio of any given type of fractures may be estimated for the given method,
2. case controlled studies, in which fracture cases and age-and-sex-matched controls are compared with respect to estimates to skeletal status obtained by the method according to the invention and optionally further to DXA-BMC measurements, or
3. by relating bio-mechanical strength of bone specimens to estimates of the skeletal status according to the invention and optionally additionally DXA-BMC-measurements of the bone in question.
Using either distal radius, fibula, tibia, ulna, a metacarpal, a phalanges, or any bone, preferably a long or a tubular bone, in which two-dimensional image data comprising information relating to the trabecular structure of the bone can be obtained, coherent values for determined bone strength, DXA-BMC and the estimate according to the present invention of the skeletal status of the vertebrate may be obtained on the basis of a calibration using eg human post-mortem bone samples or samples from other vertebrates.
The calibration of the method of the present invention may also be obtained along the above-illustrated methods as these methods generate skeletal status information and estimation on the basis of corresponding image data.
Bone strength may, e.g., be evaluated by directly measuring the plasticity and maximum load of a bone in a stress-strain diagram. (See, e.g., Lis Mosekilde). These parameters are preferably obtained in both the direction of the vertical trabeculae and in the orthogonal direction in order to have the most complete estimate of the skeletal status. However, the presently preferred manner of evaluating bone strength is using a setup determining the fracture load of the proximal femur and being designed to simulate a fall to the side of the individual.
In the method of the invention, substantially all types of models may be used for relating the chosen end point variables to a set of explanatory variables. At present, a preferred model is selected from the group consisting of: a General Linear Model, a Generalized Linear Model, an Artificial Neural Network, a Causal Probabilistic Net or Classification And Regression Trees.
As described above, the present method may be used to determine aspects of the skeletal status of any vertebrate. Thus, the vertebrate may be a human, a horse, a great ape, a large ape, an anthropoid ape, a pig, a cow, etc, and the actual bone is preferably a tubular bone, such as a bone chosen from the group consisting of radius, ulna, tibia, fibula, metacarpal, phalanges, and femur. Especially bones comprising a tube (diaphysis) and at least at one end an epiphysis comprising trabecular bone are suitable for use in the present method.